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In this paper, the honeycomb Hubbard model in optical lattices is investigated using 0(3) non- 
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and find that the narrow region could be widened by hole doping. 
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I. INTRODUCTION 

Recently, using ultracold atoms to form Bose-Einstein 
Condensates (BEC) or Fermi degenerate gases for precise 
measurements and simulations of quantum many-body 
systems, is quite impressive and has become a rapidly- 
developing field^. Experimental realizations of quan- 
tum many-body systems in optical lattices have made a 
chance to simulate strongly correlated systems^. People 
have successfully observed the Mott insulator-superfiuid 
transition both in bosonioa and fermionic atoms (e.g. 
6 Li or 40 K, etc.)^. In particular, the two-dimensional 
(fermionic) Hubbard model is one of the most interest- 
ing issues depicting the nature of high-temperature su- 
perconductivity. It is still a big challenge to clearly un- 
derstand the physics of repulsive Hubbard model on two 
dimensional (2D) lattices. Thus people try to simulate 
the Fermi-Hubbard model using a two-component mix- 
ture of repulsively interacting fermions&i'^. 

In this paper, we focus on the two dimensional Hub- 
bard model in honeycomb lattices (2D honeycomb Hub- 
bard model). In Refsi 1 ^ 11 ' 12 , it is proposed that the 
2D honeycomb optical lattice can be realized in the cold 
atoms with three detuned standing-wave lasers, of which 
the optical potential is given by 

V(x,y)= sin 2 [fc L (.t cos 6»j+y sin ^0 + ^/2] (1) 

3=1,2,3 

where 6\ = 7r/3, 82 = 27r/3, #3 = 0, and fcj, is the op- 
tical wave vector. When two-component fcrmions with 
repulsive interaction are put into the 2D honeycomb op- 
tical lattice, one can get an effective honeycomb Hubbard 
model. It is easy to change the potential barrier V by 
varying the laser intensities to tune the Hamiltonian pa- 
rameters including the hopping strength t (t-term) and 
the particle interaction U ([/-term). This lays the foun- 
dation for our later discussions and calculations. 

The simplest approach to study the honeycomb Hub- 
bard model is the Hartree-Fock (HF) mean field method, 
from which people obtain a semi-metal-insulator (MI) 
transition at a critical value (JJ/t) c between a semi-metal 
(SM) and an antiferromagnetic (AF) insulator. In the 
weak interaction region U/t < (U/t) c , the ground state 
is a semi-metal (SM) with nodal fermi-points. In the 
strong interaction region U/t > (U/t) c , the ground state 
becomes an insulator with massive fermionic excitations. 
However, the HF theory does not keep the spin rotation 




FIG. 1: Illustration of a honeycomb lattice with A and B 
sublattices. a is the length of the hexagon side and chosen to 
be unit. 



symmetry by fixing the spins along z-axis. So the re- 
sults are not reliable. Though the semi-metal-insulator 
transition of the honeycomb Hubbard model has been 
studied by different approaches ) 13 ' 14 ^ 5 ' 16 the results are 
not consistent with each other. 

In this paper, we will investigate the two dimensional 
honeycomb Hubbard model by an approach proposed in 
Refs i 17 i 18 ' 19 i 20 ' 21 i 22 that keeps spin rotation symmetry. 
By it, we find anomalous spin dynamics not far from the 
critical point of MI transition : there may exist a nar- 
row non-magnetic insulator. The narrow non-magnetic 
insulator will be in favor of hole-doping. 

The paper is organized as follows. In Sec. II, the semi- 
metal-insulator transition is studied by HF mean field 
approach. In Sec. Ill, an effective 0(3) non-linear a 
model (NLerM) is obtained to investigate properties of 
the honcycomb-Hubbard model. In Sec. IV, a global 
phase diagram is given and magnetic properties of the 
insulator state is studied based on the NLcM. In Sev. V, 
we discuss how to observe the anomalous spin dynamics 
in optical lattice of cold atoms. In Sec. VI, we show the 
doping effect on the magnetic properties of the ground 
state. Finally, the conclusions are given in Sec. VII. 
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II. METAL-INSULATOR TRANSITION 

As a starting point, the Hamiltonian of the Hubbard 
model on honeycomb lattice is given by 



h. 



Tiifriii—fj, 



E £ i £ - ( 2 ) 



Here Cj = (ci\,Ci\) and cJ CT , Cj j(T are electronic creation 
and annihilation operators, t is the hopping integral. 
U is the on-site Coulomb repulsion, a are the spin- 
indices representing spin-up(er =j) and spin-down(<r =J.) 
for electrons, /i is the chemical potential which is y in 
half-filling. denotes two sites on a nearest-neighbor 

link. rii| and are the number operators of electrons 
with up-spin and down spin respectively. 

Because the honeycomb lattice is a bipartite lattice 
(See Fig.l), we divide the system into two sublattices, A 
and B. Using the Fourier transformations, the electronic 
annihilation operators on two sublattices are written into 



Ci£A,cr 



\ " -ik Rj - 



(3) 
(4) 



where 7V S denoting the number of unit cells. For free 
fermions, the Hamiltonian could be transformed in the 
momentum space as 



k,o- 



o 6c 



flkcr 
&k<T 



where the energy of free fermions is 



ik-S 



(5) 



(6) 



U 3 + 2 cos ( V3fcy ) + 4 cos (3fc x /2) cos ( V3ky/2 
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Here the nearest neighbors of a electron in A lattice are 
defined as 

S x = | (l, V5) , 5 2 = | (l, -V3) , <y 3 = (-a, 0) (7) 

where a is the length of the hexagon side and chosen to 
be unit. The spectrum for free electrons is then obtained 
as E k = ± |£ k | • 

Next we use the path-integral formulation of electrons 
with spin rotation symmetry to study the on-site repul- 
sive interaction in Hubbard model i 17 ' 18 ' 19 ' 20 ' 21 ' 22 The in- 
teraction term can be handled by using a SU(2) invariant 
Hubbard-Stratonovich (HS) decomposition in the arbi- 
trary on-site unit vector fi^ 



(8) 



where a — (o~ x , a y ,cr z ) is the Pauli matrix. Then the HS 
transformation for the interaction term is 



i i 
iA c c\ci - Aic\toi-<JCi]). 



(9) 



Here A c and Aj are the auxiliary fields. By replacing 
electronic operators c\ a and c J;(T to Grassmann variables 
c* a and Cj jCT , the effective Lagrangian in terms of Grass- 
mann variables c* and a a is then obtained as 



C cS = E r '..- i) - r ' - ~ E' : '<•••'''< + h - c ^ _ E -^■''' l il --T<-, 

1 



E 



A2 



f7 



(A 2 + A?) + (iA c - /x) c* Cj 



(10) 



The ground state of the honeycomb Hubbard model is 
known to be long-range AF order in the large U limit. 
Such an order can be described by a simple saddle-point 
Lagrangian by fixing the direction vector field i~2^ to z- 
axis JT; = ( — l) J z and choosing the amplitude as 



iA c =-(c| Ci } = ^ 
U I t 



(11) 

(-1)*^ = (-!)* A (12) 



where n is the average on-site electron density and can- 
cels with the chemical potential [i in half-filled (n = 1) 
case. M is the staggered magnetization and A = 
is the energy band gap. The effective Lagrangian turns 
into 



£ e ff = E C *^ d rCi,a-'^2(ti,jC*C j +h.C.)-'^2 (-if Ac*(7 z 



One may obtain the spectrum of the electrons as 



£ k = ±VI6 



A 2 



(13) 



(14) 



Finally we derive the self-consistency equation for M by 
minimizing the free energy at temperature T in the Bril- 
louin zone as 



^E^ t; 



N^2E k 

k 



(15) 



Here 8 — t-^p- N is the number of the sites. 

From Fig. 2, one could find that MI transition oc- 
curs at a critical value about U/t ~ 2.23 at zero 
temperatur o 14 ' 15 i 16 i 23 . In the weakly coupling limit 
(U/t < 2.23), the ground state is a semi-metal (SM) 
with nodal fermi-points. In the strong coupling region 
(U/t > 2.23), due to M / 0, the ground state becomes 
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FIG. 2: Staggered magnetization M of half-filled honeycomb 
lattice Hubbard model at T = 0. 



an insulator with massive fcrmionic excitations. How- 
ever, the non-zero value of M only means the existence 
of effective spin moments rather than a long range AF 
order since this result is obtained in a mean field level by 
fixing the spins along z-axis. Thus one needs to examine 
stability of magnetic order against quantum fluctuations 
of effective spin moments by keeping spin rotation sym- 
metry. 



III. EFFECTIVE NONLINEAR a MODEL IN 
THE INSULATOR STATE 



In this section, we will derive an effective NLcM of spin 
fluctuations with spin rotation symmetry in the honey- 
comb Hubbard model beyond above mean field theory. 

To describe the spin fluctuations, we use the Haldane's 
mapping: 



a = (-l)Wl-L 



(16) 



where rij is the Neel vector that corresponds to the long- 
wavelength part of fti with a restriction n 2 = 1. Lj 
is the transverse canting field that corresponds to the 

short-wavelength parts of f2j with a restriction Lj ■ rij = 
0±&2!i25. We then rotate n . t0 

z-axis for the spin indexes 
of the electrons at i-site>£^&2&2L22. 

tpi = Ujci 
Ujiii ■ alii = a z 

U}U ■aU i =l i -a (17) 
where Ui sSU(2)/U(l). One then can derive the follow- 



ing effective Lagrangian after such spin transformation: 

£ e ff = ^2 ll>t d ri>i + ^2 ^* a ° W ^ 
i i 



<ij> 



(-i)vw i-n+ii- o 



*l>i (18) 



where the auxiliary gauge fields a-y = a,ij t i<j x + aij^y 
and clq (i) = ao,i (i) o x + ao,2 (i) o y are defined: 



( to j=u t Uj: a (i) = U}d T U i , 



(19) 



In terms of the mean field result M = (—1)* (ip*a z tpi) 
as well as the approximations, 



I 2 

1-1? -1-4, e % 



1 + idi 



2 i ' --hj ' 

we obtain the effective Hamiltonian as: 

CeK ^ ^ ^i d TTpi + ^ ~~ ^ ' 

i i 

- y^ihj'^ii 1 + ia-ij^j + h.c] 

m 

-A^(-l)V;^ + AM^|. 



(20) 



By integrating out the fermion fields ip* and ipi, the ef- 
fective action with the quadric terms of [arj (i) — Aer • lj] 
and ciij becomes 

1 f 13 2A 2 
5 cff = - / drJ2[~^(ao (i)-Aa ■ h) 2 +Ap s a 2 ij +—l 2 i ) 

(21) 

where p s and <; are two parameters. 

To learn the properties of the low energy physics, 
we study the continuum theory of the effective action 
in Eq. (|21[) . In the continuum limit, we denote rij, lj, 
idij ~ UjUj - 1 and a (£) = U}d T Ui by n(x,y), l(x,y), 
Wd x U (or WdyU) and Wd r U, respectively. Fr om the 
relations between Wd^U and 9 M n, 



2 

a T = 


2 

a r,l " 


2 

- a r,2 = 


< = 


a l,l ~ 


^ <2 = 


ao-l=- 


-— (n x <9 T n 

2 v 



(22) 



the continuum formulation of the action in Eq. (|21[) turns 
into 



(d M n) 2 , M = ^2/, (23) 
(24) 



ijf dr J d 2 r[<;(d T n) 2 + ft ( V nf 



4zA<j (n x <9 T n) • 1 + ( 



2A 2 
17 



4A 2 01 2 



(25) 
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where the vector ao is defined as ao = (ao,i, Q.0.2, 0) . 

Finally we integrate the transverse canting field 1 and 
obtain the effective NLcrM as 



S cB = ^ I dr f d 2 r[- (d T nf + c ( V n) 2 ] 
25 Jo J c 



(26) 



with a constraint n 2 = 1. The coupling constant g and 
spin wave velocity c arc defined as: 



c 2 Ps 
9=—,c = — , 

Ps X 
Here p s is the spin stiffness, 



1 



N V4(|e k | 2 + A2)l 

where the corresponding coefficient e 2 is 
1 



(27) 



(28) 



4 



27r 



2A 2 + 27i 2 ) cos ( VZk, 



+ 36i 2 cos (3k x /2) cos ^v / 3fc J) /2j cos y\/3k y 
+ 2 (5A 2 + 27t 2 ) cos (3/^/2) cos (V3ky/2 



9t 2 cos(3k x ) (l + cos (Vzky))} 



(29) 



X 1 ' is the transverse spin susceptibility, of which ( is 



1 E 



N - 

k 4 |Ck| 2 + A 2 



(30) 



One can see detailed calculations of p s and £ in the ap- 
pendix. 

The numerical results of p s and c are illustrated in 
Fig.3, where one can find that c = 0.168901985i = 
1.055637 J in the strongly coupling limit (U = 25t) match 

the earlier results, c = 1.06066 J (J = ^-), obtained from 
Heisenberg model 2 ^. 

In addition, we need to determine another important 
parameter - the cutoff A. On the one hand, the effec- 
tive NLcrM is valid within the energy scale of Mott gap, 
2A = UM. On the other hand, the lattice constant is a 
natural cutoff. Thus the cutoff is defined as the following 
equation^ 

2A 

A = min(l, ) (31) 

c 



IV. MAGNETIC PROPERTIES OF THE 
INSULATOR STATE 

In this section we will use the effective NLcrM to study 
the magnetic properties of the insulator stated. The 
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FIG. 3: The spin stiffness p 3 and the spin wave velocity c of 
the half- filled honeycomb Hubbard model at T = 0. 




FIG. 4: The dimensionless coupling constant a of Hubbard 
model on honeycomb lattice of half filling at T = 0. Insert 
(b) shows non-magnetic region. 



Lagrangian of NLcrM with a constraint (n 2 = 1) by a 
Lagrange multiplier A becomes 



C 



off 



1 

2^9 



(cLn) 2 + C 2 ( V n) 2 +iA(l-n 2 ) (32) 



where iX = m 2 and m is the mass gap of the spin fluctu- 
ations. 

Using the large- N approximation we rescale the field 
n -^\fNn and obtain the saddle-point equation of motion 
as 

(n ) 2 -k B T ]T n(q,iw n ) = l. (33) 

In Eq. (|33|) . no is the mean field value of n and 
II(q, itu n ) = - ji +c Q i+m2 is the propagator of the spin 
fluctuations Sn = n— no. Here lo„ = 2imk^T, n = inte- 
gers. 

At finite temperature, the solution of no is always zero 
that is consistent to the Mcrmin-Wigncr theorem. From 
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FIG. 5: The mass gap m of the spin fluctuations of Hubbard 
model on honeycomb lattice of half filling at T = 0. 



Eq. (f33|) . we may get the solution of m as 

cA 



e gk B T smb. 



2k B T 



In the limit T <C A, Eq. (|34|) can be rewrite as: 



m = 2fcB7 1 sinh ' ■{ — e? 



2ttc /l 1 



^bT \g g c 



where 



9c 



4ir 

T' 



(34) 



(35) 



(36) 



Therefore, at zero temperature the solutions of no and m 
of Eg . (|3"3"]) are determined by the dimcnsionless coupling 
constant a = gA. In particular, there exists a critical 
point a c = 4n (or g c = ^): For the case of a < 47r, we 
get solutions of rig and m: 



n = (l- A)i/2 jm = 
3c 



(37) 



For the case of a > 47r, we get solutions of no and m: 

1 1. 



np = 0, m = 47rc(- 



5c g 



(38) 



So we calculate the dimcnsionless coupling constant 
a = gA and show results in Fig. 4. The quantum critical 
points corresponding to a c ~ 4tt turn into (U/t) c2 ~ 2.88 
and (U/t) c3 ~ 2.93 which divide the insulator state into 
three phases - a quantum disordered state (QD) in the 
region of 2.88 < U/t < 2.93 and two long range AF order 
in the regions of 2.23 < U/t < 2.88 and U/t > 2.93. 

In the regions of 2.23 < U/t < 2.88 and U/t > 2.93 
(where a < a c ), at low temperature the mass gap m of 
spin fluctuations is determined by: 



2fcBTexp 



2ttc (\ 



1 



(39) 



FIG. 6: Spin order parameter Mo of the Hubbard model on 
honeycomb lattice of half filling at T — 0. Insert (b) shows 
non-magnetic region. 



Because the energy scale of the mass gap m is al- 
ways much smaller than the temperature, i.e., m <C 
k^T (or uj n ), quantum fluctuations become negligible 
in a sufficiently long wavelength and low energy regime 
(m < |cq| < ksT) . Thus in this region one may only con- 
sider the purely static (semiclassical) fluctuations. 

At zero temperature, the mass gap m vanishes (See 
Fig. 5), which means that long range AF order appears. 
To describe the AF order, we introduce a spin order 



parameter* 



28,29.30. 



Mr 



M M , 0,1/0 , 
2 2 g c 



As shown in Fig. 6, the ground state of AF ordered phase 
has a finite spin order parameter. 

In the region of 2.88 < U/t < 2.93 (where a > a c ), 
there is a finite mass gap of spin fluctuations at zero 
temperature (See Fig. 5): 



4ttc( 



1 



i) 

9 



(41) 



Therefore, the ground state of the insulator in this region 
is not a long range AF order. Instead, it is a quantum 
disordered state (or non-magnetic insulator state) with 
zero spin order parameter M = in a narrow non- 
magnetic window (See Fig. 6). 

Based on above results, we get the global phase dia- 
gram which is illustrated in Fig. 7. One can see that at 
finite T, there are four crossover lines, Thf, T p , T u , that 
separate five regions. 

The highest crossover line is Thf that is obtained from 
Eq. (fT5"l) and denotes the establish of the effective spin- 
moments. Above Thf, it is metal phase without en- 
ergy gap A = 0. The crossover line T p ~ p s denotes 
the validity of the NLerM, where p s is the energy scale 
of spin stiffness. In the region T p < T < Thf, the 
free spin-moments are established (denoted by M ^ 0) 
that show a Curie- Weiss behavior. In this region one 
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FIG. 7: Phase diagram of the Hubbard model on honeycomb 
lattice of half filling at finite temperature. Insert (b) shows 
non-magnetic region. 



cannot use the effective NLtrM. Below T p , short range 
spin-correlation exists and the effective NLtrM is valid. 
The region below T p is dominated by the crossover lines 
T v ~ p s \l — g/g c \ together with two quantum critical 
points (QCP) at (U/t) c2 = 2.88 and (U/t) c3 = 2.93 which 
represent the QCP of g = g c in the NLtrM.^22 



V. SPIN-SPIN CORRELATIONS 

In order to make our theoretical predictions ver- 
ifiable, we discuss the detection of the anoma- 
lous spin dynamics in ultracold atom experiments 
via spatial spin-spin correlations (S (r, t) ■ S (0, 0)) = 
(ct (r, t) trc (r, t) ■ c 1 (0, 0) trc (0, 0)) or dynamic spin sus- 
ceptibility x" (q, w )- Here x" (^i^) i s defined as 

X " (q, w) = \ (1 - e-"?) J dtdve^-^ (S (r, t) ■ S (0, 0)) . 

(42) 

On the one hand, people may observe spatial spin- 
spin correlations (S z (ri) S z (r 2 )) by noise correlation 
spectroscop y 31 ' 32 . As demonstrated in Refi 3 -^, using a 
probe laser beam which goes through the system, one can 
measure the phase shift or change of polarization of the 
outgoing beam (A ou t) oc (M z ) to obtain the magnetiza- 
tion (M z ) oc J dvcj){r) (S z (r)), where (f>(r) is the spatial 
intensity profile of the laser beam. The quantum noise 
(l D 2 ut ) oc {Ml) a/dndra^CrO^raX^CrOX^Cra)) 
reveals the atomic correlations in the system. This 
method can directly examine the existence of 2D AF cor- 
relations. 

On the other hand, people may also observe dynamic 
spin susceptibility x" (l) w ) which links with spin-spin 
correlations by NMR measurement. Consequently, it de- 
termines the spin-spin correlation in the real space as 
cos (Qo • r) e~l r l/'>, with the correlation length £ and the 
AF wave- vector Qo = (±-,±— ). To show clearly the 



FIG. 8: Correlation lengh £ at T = 0.02, 0.1, 0.2. The unit 
of temperature is t. 



anomalous spin dynamics, we calculate the spin-spin cor- 
relation (S (r, t) ■ S (0, 0)) . The dynamic spin susceptibil- 
ity at Qo is determined by n fields 

gc 

X(q,w) = (n a (q,o;)n (-q,-u/)} = -r- — 5- ^ (43) 

where we have introduced the momentum q. The 
equal-time spin-spin correlation function is proportion 
to e~' r '/£ where the spin-correlation length £ is there- 
fore given by — . As Eq. (|4"2")) indicates, the spin dynamic 
structure factor as well as the dynamic spin susceptibility 
function may reflect the effective short-range magnetic 
correlation length £. 

In Fig. 8, we give the spin-correlation length defined as 
£ = -L at T = 0.02i, O.lt, 0.2t, respectively. From the 
Fig. 8, taking T = Q.02t as an example, one can see that 
the spin-correlation length increases quickly with increas- 
ing interaction U/t. However, the spin-correlation length 
doesn't increase monotonously with U/t - it will de- 
crease and reach a minimum value near the MI transition 
U/t ~ (U/t) c2 ~ (U /t) c3 . The dip of the spin-correlation 
length will indicate the existence of the non-magnetic 
state near MI transition. When one increases the inter- 
action further, the spin-correlation increases again and 
finally decreases in the strongly interacting limit due to 
J — > 0. For other cases with higher temperature, T = 
O.lt, 0.2i, there exist similar dip structure of the spin- 
correlation length via U/t which means that people may 
observe the anomalous spin dynamics more easily in ex- 
periments. 



VI. DOPING EFFECT 

In this section we will leave from half-filling and study 
the hole doping cases. In case of hole-concentration d = 
1 — n, the chemical potential [i is not -j any more. The 
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FIG. 9: Spin stiffness p s of the Hubbard model on honeycomb 
lattice with hole-concentration d = 0.0, 0.1, 0.2, 0.3, 0.4, 0.5 
at T = 0. 



FIG. 10: Spin wave velocity c of the honeycomb lattice Hub- 
bard model with hole-concentration d = 0.0, 0.1, 0.2, 0.3, 0.4, 
0.5 at T = 0. 



Hamiltonian becomes 



H = - •:/,., '-,'(% + h.c.) - ^2 AjC*fV<7Cj 



(ij) 



E 



Un 



(44) 



At this case, we could obtain M, p, £, p s similarly with 
that of half-filling case in T = as follows, 



1 = 



A 2 



L E 



AT 



^<^4(iar + A ; 

c 2 



7V ^M 4 (I^| 2+A2 ) 



Given certain hole-concentration d, the p s and c may 
be obtained. From comparation in Fig. 9, one can see 
that spin stiffness p s rises when the hole-concentration 
increases from d = 0.1 to d = 0.5, and it shows a dif- 
ferent behavior in the large doping density from that of 
low doping density. From Fig. 10, the spin wave velocity 
c goes up as well when we increase hole concentration 
and it becomes much larger at critical points. The di- 
mcnsionless coupling constant a are calculated in Fig. 11 
in different hole-concentrations, from which one can see 
that all critical points go up as we increase doping con- 
centration d. And more importantly, we find that the 
non-magnetic regions widen when the hole-concentration 
increases. 

In order to show the doping effect clearly, we plot the 
phase diagram of the honeycomb lattice Hubbard model 
at T = (See Fig. 12). From the phase diagram, one can 
see that there are three critical lines (U c i, U c i-, U c s) sep- 
arating four regimes - one SM regime, two AF regimes 



a 




FIG. 11: The dimensionless coupling constant a of the Hub- 
bard model on honeycomb lattice with hole-concentration d = 
0.0, 0.1, 0.2, 0.3, 0.4, 0.5 at T = 0. 



and one non-magnetic regime. In Fig. 12, U c \ that has 
been given in Ref.— is derived by Hartree-Fock approxi- 
mation. U C 2 and U C 3 are determined by the critical point 
of the effective NLtrM. In particular, the region of quan- 
tum disordered state becomes much wider by increasing 
the hole-concentration. 



VII. CONCLUSION 

In this paper, we investigate the two dimensional hon- 
eycomb Hubbard model with an approach that keeps spin 
rotation symmetry. By it, we found anomalous spin dy- 
namics not far from the critical point of MI transition : 
there may exist a narrow non-magnetic insulator. The 
narrow non-magnetic insulator will be in favor of hole- 
doping. Also we discuss how to observe the anomalous 
spin dynamics in cold atoms. 

In the end we discuss the nature of the non-magnetic 
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FIG. 12: Phase diagram of the Hubbard model on honey- 
comb lattice at different hole-concentration at T — 0. The 
square dot line, circle dot line, triangle dot line represent 
Ud, Uc2, Uc3, respectively. 



insulator in the two dimensional honeycomb Hubbard 
model. In Refi^, it is pointed out that by considering the 
fermionic nature of vortices (half-skyrmions), nodal spin 
liquid becomes the ground state of the non-magnetic in- 
sulator state in the 7r-flux Hubbard model on square lat- 
tice. There exist three types of quasi-particles in nodal 
spin liquids: nodal fermionic spinons, gapped bosonic 
spinons and roton-like U(l) gauge field. Following the 
similar approach, one may draw the same conclusion that 
the non-magnetic insulator in the two dimensional honey- 
comb Hubbard model is another example nodal spin liq- 
uid with similar quasi-particles - nodal fermionic spinons, 
gapped bosonic spinons and roton-like U(l) gauge field. 
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APPENDIX A: THE DETAILED CALCULATIONS 

OF p s AND ? 

To give p s and q for calculation, we choose U% to be 

(Al) 



u t = ( f 



where = Zjcrz,, Zj = (zif, Zjj.) , z^z^ 1. And the spin 
fluctuations around n,- = z,- is 



rij = Zj+Re (fa) x+Im (&) y 
l-|^| 2 /8^ +0 ^ 



(A2) 
(A3) 



Then the quantities UjlJj and Ujd T Ui can be expanded 
in the power of fa — <f>j and d T (j>i, 



UjUj = e^^^y 

tj\qtj.-( , 

t/ '^'" o 



(A4) 
(A5) 



According to Eq. (fT9|l . the gauge field ay and ao (i) are 
given as 

1 

2 



«o (i) = -^d T 4>i<T y . 



(A6) 
(A7) 



Supposing dij and ao (i) to be a constant in space and 
denoting di<f>i= a and d T (f>i = iB y , we have 



(A8) 
(A9) 



"0 (*) = -^ByCTy. 

The energy of Hamiltonian of Eq. (f2"Tj) becomes 

E(B y ^) = \^Bl + ^ PsS ?. (A10) 

Then one could get £ and p s from the following equations 
by calculating the partial derivative of the energy 



c 



1 d 2 E Q (By ) 

N dBl 
1 d 2 E (a) 



\By=0 



s N da 2 



(All) 



(A12) 



Here Eq (B y ) and Eq (a) are the energy of the lower Hub- 
bard band 



^(*v) = £(<k + <k) (A13) 

k 

E ( a )=J2(KM + E -M) ( A14 ) 

k 

where £7^. k , E_ k and £^ k , E p _ k are the energies of the 
following Hamiltonian "w- and Ti p 

<»J> i 

+ £V*ao(<)^ (A15) 



<jj> 



(A16) 

Using the Fourier transformations for W» , we have the 
spectrum of the lower band of H.^ 



4, k = -\/fl&l±^) 
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where £k has been obtained in Eq.©. Using Eq. (|All[ ). Using Eq. (|A12p . p s is given as 
£ is obtained as 

-r. 



N ^ ( 

k 4(iai 2 + A2 



(A18) 



Ps = 



N V4(|e k | 2 + A2)4 



(A22) 



The corresponding coefficient e 2 is 



Similarly, using the Fourier transformations for 7i p , we 
obtain the spectrum of the lower band of H p 



El 



±,k 



A 2 + \ip\ 2 + i^r ± 4A 2 - - i><p*Y 



where ip and ip arc defined as 

<p = -tY^e lkS cos Qa- 5 

= -i^e ik " 5 sin (-a- 5 



(A19) 

(A20) 
(A21) 



e 2 = -t 2 [6A 2 + 27t 2 + (2 A 2 + 27i 2 ) cos 



+ 36f 2 cos (ik x /2) cos {Viky/2^ cos (Viky^j 
+ 2 (5A 2 + 27£ 2 ) cos (3fc x /2) cos (V3k y /2\ 
+ 9i 2 cos (3^) (l + cos (V3ky^}. 



(A23) 
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